/*
  Author: Kuljit Virk
  Date  : 22 Jan 2012
  Purpose: Input output for the superlattice code

  Read in: Parameter file for tight-binding Hamiltonians
*/
#include "tbparameters.hpp"

using namespace arma;
using namespace std;


void hamiltonian(struct tb p,rowvec kvec,cx_mat & H,cx_mat & Hc,cx_mat & Ha,
		 cx_mat & Vup,cx_mat & Vdn){
  
  float t=1.0;
  cx_rowvec edk(4),gdn(4),gup(4);
  cx_mat HR,HL;
  mat d;
  
  // d1    d2    d3   d4
  // 1     1    -1    -1;  x
  // 1    -1     1    -1;  y 
  // 1    -1    -1     1;  z

  d << 1.0  <<   1.0  <<  -1.0  <<  -1.0 << endr
    << 1.0  <<  -1.0  <<   1.0  <<  -1.0 << endr
    << 1.0  <<  -1.0  <<  -1.0  <<   1.0 << endr;

  edk  = pi/2 * cx_double(0.0,1.0) * (kvec * d);
  edk  = exp(edk);
    
  gdn << edk(1) + edk(2) << edk(1) - edk(2) << -edk(1) + edk(2) << -edk(1) - edk(2) << endr;
  gup << edk(0) + edk(3) << edk(0) - edk(3) <<  edk(0) - edk(3) <<  edk(0) + edk(3) << endr;
  
  gup = gup/4;
  gdn = gdn/4;

//   /*% up means Anion plane above Cation in physical space
//     % In Hamiltonian matrix, it is the row below. 
//     % Output of Ham is [Ha Vac;Vac' Hc] where a = anion and c = cation
//     %
//     % So for up: [Ha Vac;Vac' Hc] --> [Hc Vac';Vac Ha]
//   */
  Ham(p,gup,H, HL, HR, Vup);
  Vup = Vup.t();
  Hc  = HR;
  Ha  = HL;
    
  // So for dn: [Ha Vac;Vac' Hc] remains as is
  cx_mat Htemp;
  Ham(p,gdn,Htemp, HL, HR, Vdn);

  int n=Vup.n_rows;
  
  H.submat(0,0,n-1,n-1)     = Hc;
  H.submat(n,n,2*n-1,2*n-1) = Ha;

  H.submat(0,n,n-1,  2*n-1) = Vup+Vdn.t();
  H.submat(n,0,2*n-1,  n-1) = Vup.t()+Vdn;

  return;
};

void velocitymatrix(struct tb p,rowvec kvec,cx_mat & H,cx_mat & Hc,cx_mat & Ha,
		    cx_mat & Vup,cx_mat & Vdn,int cartesiancoord){
  
  float t=1.0;
  cx_rowvec edk(4),gdn(4),gup(4);
  cx_mat HR,HL;
  mat d;
  int i=cartesiancoord;
  // d1    d2    d3   d4
  // 1     1    -1    -1;  x
  // 1    -1     1    -1;  y 
  // 1    -1    -1     1;  z

  d << 1.0  <<   1.0  <<  -1.0  <<  -1.0 << endr
    << 1.0  <<  -1.0  <<   1.0  <<  -1.0 << endr
    << 1.0  <<  -1.0  <<  -1.0  <<   1.0 << endr;

  edk  = pi/2 * cx_double(0.0,1.0) * (kvec * d);
  edk  = exp(edk);

  // Replace each exp(ikd) by its derivative id(n) exp(ikd)
  // where n = cartesian coordinate. Note that k = kvec * pi/2 so
  // the differentiation is with respect to (kvec*pi/2) so that only
  // d is produced. Remember to scale the velocity matrix by the scale
  // of length, which is the lattice spacing, in the final result.
  for (int j=0;j<4;j++)
      edk(j) = cx_double(0.0,1.0)*d(i,j)*edk(j);
    
  gdn << edk(1) + edk(2) << edk(1) - edk(2) << -edk(1) + edk(2) << -edk(1) - edk(2) << endr;
  gup << edk(0) + edk(3) << edk(0) - edk(3) <<  edk(0) - edk(3) <<  edk(0) + edk(3) << endr;
  
  gup = gup/4;
  gdn = gdn/4;
  
//   /*% up means Anion plane above Cation in physical space
//     % In Hamiltonian matrix, it is the row below. 
//     % Output of Ham is [Ha Vac;Vac' Hc] where a = anion and c = cation
//     %
//     % So for up: [Ha Vac;Vac' Hc] --> [Hc Vac';Vac Ha]
//   */
  Ham(p,gup,H, HL, HR, Vup);
  Vup = Vup.t();
  Hc  = HR;
  Ha  = HL;
  
  // So for dn: [Ha Vac;Vac' Hc] remains as is
  cx_mat Htemp;
  Ham(p,gdn,Htemp, HL, HR, Vdn);
  
  int n=Vup.n_rows;
  
  H.zeros();

  // See Table A in Vogl and observe that only terms with k-dependent
  // factors are the cation-anion off-diagonal terms
  H.submat(0,n,n-1,  2*n-1) = Vup+Vdn.t();
  H.submat(n,0,2*n-1,  n-1) = Vup.t()+Vdn;

  // Set the diagonal to zero because it does not have any
  // k dependent terms and the g(k) factors are not included
  // when storing the diagonal elements. So must do this manually
  H.diag(0)  *= 0.0;
  Hc.diag(0) *= 0.0;
  Ha.diag(0) *= 0.0;

  return;
};



void Ham(struct tb p,const cx_rowvec& g1,
	 cx_mat& H, cx_mat& HL, cx_mat& HR, cx_mat& VLR){

  bool sstar;
  cx_rowvec g(5);
  g(1)=g1(0);
  g(2)=g1(1);
  g(3)=g1(2);
  g(4)=g1(3);

  // Without SO coupling, single spin matrix
  // H diagonal: Increase dim by 1 to use 1-indexing
  // for direct conversion from matlab. discard row 0 and col 0 at end
  cx_vec E;
  E << 0
    << p.Esa
    << p.Epa
    << p.Epa
    << p.Epa
    << p.Esc
    << p.Epc
    << p.Epc
    << p.Epc;

  H = diagmat(E);
          
  // H(sa,sc)
  H(1,5) = p.Vss*g(1);
  // H(sa,pc)
  H(1,6) =  p.Vsapc*g(2);
  H(1,7) =  p.Vsapc*g(3);
  H(1,8) =  p.Vsapc*g(4);
  // H(pa,sc)
  H(2,5) = -p.Vscpa*g(2);
  H(3,5) = -p.Vscpa*g(3);
  H(4,5) = -p.Vscpa*g(4);
        
  // H(pa,pc) same orbital
  H(2,6) = p.Vxx*g(1);
  H(3,7) = p.Vxx*g(1);
  H(4,8) = p.Vxx*g(1);
        
  // H(pa,pc) cross orbitals
  H(2,7) = p.Vxy*g(4);
  H(2,8) = p.Vxy*g(3);
  H(3,6) = p.Vxy*g(4);
  H(3,8) = p.Vxy*g(2);
  H(4,6) = p.Vxy*g(3);
  H(4,7) = p.Vxy*g(2);
                
  cx_mat h1=H.submat(1,1,8,8);

  if (strcmpi(p.name,"vogl") || strcmpi(p.name,"klimeck")){
    sstar=true;
    H.zeros(11,11);
    H.submat(1,1,8,8)=h1;  
    H(9,9)   = p.Essa;
    H(10,10) = p.Essc;

    for (int i=2;i<=4;i++){
      H(i,10)  = -p.Vpassc*g(i);
      H(i+4,9) =  p.Vssapc*g(i);
    }    
  }
  else
    sstar=false;
  
  //-----------------------------------------------------------------------
  // SO Coupling: if applicable based on p.name
  //-----------------------------------------------------------------------
  // n = index at the end of the second block
  // but since n_rows is 1 greater than the relevant rows
  // set n as:
  int n = H.n_rows-1;
  if (strcmpi(p.name,"klimeck") || strcmpi(p.name,"chadi")){
    //// Double the matrix for SO coupling
    h1 = H.submat(1,1,n,n);  
    H.zeros(2*n+1,2*n+1);
    H.submat(1,1,n,n) = h1;
    H.submat(n+1,n+1,2*n,2*n) = h1;
    //// Add SO coupling: From Harrison's book            
    // ANION
    so(H,p.soa,n,2,3,4);
    // CATION
    so(H,p.soc,n,6,7,8);
  }
  //-----------------------------------------------------------------------
  //// Symmetrization
  //-----------------------------------------------------------------------
  E           = H.diag();
  cx_mat    V = H - diagmat(E);
  H           = diagmat(E) + V + V.t();


  n = H.n_rows-1;

  //-----------------------------------------------------------------------
  // Rearrangement to put Anion and Cation orbitals together 
  //-----------------------------------------------------------------------
  // 1. Index permutation array to put the s* orbitals with anion and cation blocks.
  irowvec idx(n+1);
  if (sstar)
    if (n==10)
      idx << 1 << 2 << 3 << 4 << 9 << 5 << 6 << 7 << 8 << 10;
    else
      idx << 1 << 2 << 3 << 4 << 9 << 5 << 6 << 7 << 8 << 10 
	<< 11 << 12 << 13 << 14 << 19 << 15 << 16 << 17 << 18 << 20;
  else
    for (int i=0;i<n+1;i++)
      idx(i) = i+1;

  h1.zeros(n,n);
  h1 = H.submat(1,1,n,n);
  for (int i=0;i<n;i++)
    for (int j=0;j<n;j++){
      H(i+1,j+1)=h1(idx(i)-1,idx(j)-1);
    }
  // 2. Index permutation array to put both spins of anion and cation in adjacent blocks
  if (n==20)
    idx <<  1 <<  2 <<  3 <<  4 << 5 
	<< 11 << 12 << 13 << 14 << 15
	<<  6 <<  7 <<  8 <<  9 << 10
	<< 16 << 17 << 18 << 19 << 20;
    else
      if (n==16)
	idx <<  1 <<  2 <<  3 <<  4 
	    <<  9 << 10 << 11 << 12 
	    <<  5 <<  6 <<  7 <<  8 
	    << 13 << 14 << 15 << 16;
      else
	for (int i=1;i<n;i++)
	  idx(i-1) = i;
  // 3. Apply the index permutation to the array
  h1 = H.submat(1,1,n,n);
  for (int i=0;i<n;i++)
    for (int j=0;j<n;j++){
      H(i+1,j+1)=h1(idx(i)-1,idx(j)-1);
    }
  h1  = H;
  //-------------------------------------------------------------------------
  // Now just discard the 0 row and 0 column
  // And read off matrices the matlab way
  H   = h1.submat(1,1,n,n);
  
  HL  = h1.submat(1,1,n/2,n/2);
  HR  = h1.submat((n/2+1),(n/2+1),n,n);
  VLR = h1.submat(1,(n/2+1),n/2,n);

};

void so(cx_mat & H,
	const double lambda,const int n,
	const double xp, 
	const double yp, 
	const double zp){
  
  double xm = xp+n;
  double ym = yp+n;
  double zm = zp+n;
  
  H(xp,yp) =  cx_double(0.0,lambda);    // (px+,py+)
  //         H(xm,ym) =  cc(H(xp,yp));
  H(xm,ym) =  -H(xp,yp);        // (px-,py-)
  
  H(xp,zm) =  lambda;           // (px+,pz-)
  //         H(xm,zp) = -H(xp,zm);
  H(xm,zp) = -H(xp,zm);         // (px-,pz+)
  
  H(yp,zm) =  cx_double(0.0,lambda);    // (py+,pz-)
  //         H(ym,zp) =  -cc(H(yp,zm));
  H(ym,zp) =  H(yp,zm);         // (py-,pz+)
};
